#include <iostream>
#include <cmath>
#include <cstdio>
#include <pthread.h>
#include <stdlib.h>

#define NUM_THREADS 8

pthread_t thread[NUM_THREADS];

unsigned short Nvar1_seed[3] = {4,7,2};
unsigned short Nvar2_seed[3] = {5,8,1};
unsigned short Nvar3_seed[3] = {9,6,3};

double *trials;

int _M; // numero de iterações

pthread_mutex_t mutexsum;

long long _sum_x[NUM_THREADS], _sum_x2[NUM_THREADS];

double Normal(float mean, float std)
{
    double var_u1, var_u2, var_u3;
    double x, tmp_val;
    
    var_u1 = erand48(Nvar1_seed);
    var_u2 = erand48(Nvar2_seed);
    x = -log(var_u1);
    tmp_val = exp(-((x - 1)*(x - 1)) / 2);
    while(var_u2 > tmp_val)
    {
        var_u1 = erand48(Nvar1_seed);
        var_u2 = erand48(Nvar2_seed);
        x = -log(var_u1);
        tmp_val = exp(-((x - 1)*(x - 1)) / 2);
    }
    var_u3 = erand48(Nvar3_seed);
    if(var_u3 > 0.5)
        return(mean + (std * x));
    else
        return(mean - (std * x));
}

typedef struct
{
    float S;
    float E;
    float R;
    float SIGMA;
    float T;
    int id;
} thread_args;

double max(double a, double b)
{
    return a > b ? a : b;
}

void* blackScholes_run(void* args)
{
    thread_args *arg = (thread_args*) args;
    double t;
    
    _sum_x[arg->id] = 0;
    _sum_x2[arg->id] = 0;
    
    for(int i=0; i < _M/NUM_THREADS; i++)
    {
        t = arg->S * expl((arg->R - 0.5 * arg->SIGMA * arg->SIGMA) * arg->T + arg->SIGMA * sqrt(arg->T) * Normal(0,1));
        trials[i] = expl((-arg->R) * arg->T) *  max(t - arg->E, 0);

        _sum_x[arg->id] += trials[i];
        _sum_x2[arg->id] += trials[i] * trials[i];
    }
    
    pthread_exit((void*) arg->id);
}

int main(int argc, const char * argv[])
{
    int i, rc;
    double mean, stddev, confwidth, confmin, confmax;
    long long sum_x = 0, sum_x2 = 0;
    void *status;
    
    thread_args args;
    
    scanf("%f", &(args.S));
    scanf("%f", &(args.E));
    scanf("%f", &(args.R));
    scanf("%f", &(args.SIGMA));
    scanf("%f", &(args.T));
    scanf("%d", &_M);

    trials = (double*) malloc(_M * sizeof(double));
    
    thread_args *t_arg = (thread_args*) malloc(NUM_THREADS * sizeof(thread_args));
    
    for (i = 0; i < NUM_THREADS; i++)
    {
        t_arg[i] = args;
        t_arg[i].id = i;
        rc = pthread_create(&thread[i], NULL, blackScholes_run, (void *) &t_arg[i]);
    }
    
    for(i=0; i<NUM_THREADS; i++) {
        rc = pthread_join(thread[i], &status);
        if (rc) {
            printf("ERRO; codigo de retorno da pthread_join() : %d\n", rc);
            exit(-1);
        }
    }

    for(i=0;i< NUM_THREADS;i++)
    {
        sum_x += _sum_x[i];
        sum_x2 += _sum_x2[i];
    }
    
    mean = sum_x/_M;
    stddev = sqrt(sum_x2/_M - mean * mean);
    confwidth = 1.96f * stddev/sqrt(_M);
    confmin = mean - confwidth;
    confmax = mean + confwidth;
    
    printf("S        %.0lf\n", args.S);
    printf("E        %.0lf\n", args.E);
    printf("r        %.0lf\n", args.R);
    printf("sigma    %.0lf\n", args.SIGMA);
    printf("T        %.0lf\n", args.T);
    printf("M        %d\n", _M);
    printf("Confidence interval: (%.4lf, %.4lf)\n", confmin, confmax);
    
    free(trials);
    free(t_arg);
    
    return 0;
}

